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Abstract 

In this paper, we present results of a discontinuous Galerkin (DC) scheme applied to 
deterministic computations of the transients for the Boltzmann-Poisson system describing 
electron transport in semiconductor devices. The coUisional term models optical-phonon 
interactions which become dominant under strong energetic conditions corresponding to 
nano-scale active regions under applied bias. The proposed numerical technique is a finite 
element method using discontinuous piecewise polynomials as basis functions on unstructured 
meshes. It is applied to simulate hot electron transport in bulk silicon, in a silicon n~^-n-n~^ 
diode and in a double gated 12nm MOSFET. Additionally, the obtained results are compared 
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to those of a high order WENO scheme simulation and DSMC (Discrete Simulation Monte 
Carlo) solvers. 



Keywords: Deterministic numerical methods, Discontinuous Galerkin schemes, Boltz- 
mann Poisson systems. Statistical hot electron transport. Semiconductor nano scale devices. 
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1 Introduction 



The evolution of the electron distribution function /(t, x, k) in semiconductors in dependence 
of time t, position x and electron wave vector k is governed by the Boltzmann transport 
equation (BTE) [231 [271 [19] 

^ + ^Vke ■ Vx/ - |e ■ Vk/ = Qif) , (1.1) 

where h is the reduced Planck constant, and q denotes the positive elementary charge. The 
function e(k) is the energy of the considered crystal conduction band measured from the 
band minimum; according to the Kane dispersion relation, e is the positive root of 

. . fx hj , , 

e{l + ae) = —-^, 1.2 
2m* 

where a is the non-parabolicity factor and m* the effective electron mass. The electric field 
E is related to the doping density Nd and the electron density n, which equals the zero-order 
moment of the electron distribution function /, by the Poisson equation 

Vx [e.(x) V^V] = ^ [n{t, x) - Nd{^)] , E = -V^V , (1.3) 

where eo is the dielectric constant of the vacuum, er(x) labels the relative dielectric function 
depending on the material and V the electrostatic potential. The collision operator Q{f) 
takes into account acoustic deformation potential and optical intervalley scattering [291 [28] . 
For low electron densities, it reads 

Q(/)(t,x,k)= / [5(k',k)/(t,x,k')-5(k,k')/(t,x,k)]rfk' (1.4) 

with the scattering kernel 

^(k,k') = [ug + 1) K S{e{k') - e{k) + hup) 

+ ngK6{e{k') - ^(k) - hujp) + Ko <5(e(k') - ^(k)) (1.5) 

and K and Kq being constant for silicon. The symbol S indicates the usual Dirac distribution 
and Up is the constant phonon frequency. Moreover, 
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is the occupation number of phonons, ks is the Boltzmann constant and Ti is the constant 
lattice temperature. 

Semiclassical description of electron flow in semiconductors thus is an equation in six 
dimensions (plus time if the device is not in steady state) for a truly 3-D device, and four 
dimensions for a 1-D device. This heavy computational cost explains why the BP system 
is traditionally simulated by the Direct Simulation Monte Carlo (DSMC) methods [22j. In 
recent years, deterministic solvers to the BP system were proposed [HI EH [3], Ej HI 161120]. 
These methods provide accurate results which, in general, agree well with those obtained from 
Monte Carlo (DSMC) simulations, often at a fractional computational time. Moreover, they 
can resolve transient details for the pdf, which are difficult to compute with DSMC simulators. 
The methods proposed in [U [6] used weighted essentially non-oscillatory (WENO) finite 
difference schemes to solve the Boltzmann-Poisson system. The advantage of the WENO 
scheme is that it is relatively simple to code and very stable even on coarse meshes for 
solutions containing sharp gradient regions. A disadvantage of the WENO finite difference 
method is that it requires smooth meshes to achieve high order accuracy, hence it is not very 
fiexible for adaptive meshes. 

On the other hand, motivated by the easy /ij»-adaptivity and simple communication 
pattern of the discontinuous Galerkin (DC) methods, researchers have worked on developing 
the DC method for solving the Boltzmann equation and its macroscopic models [3, [8], |2ll 
[25l [21] . The type of DG method that we will discuss here is a class of finite element 
methods originally devised to solve hyperbolic conservation laws containing only first order 
spatial derivatives, e.g. [TU [13], [121 HB [IS]- Using completely discontinuous polynomial 
space for both the test and trial functions in the spatial variables and coupled with explicit 
and nonlinearly stable high order Runge-Kutta time discretization, the method has the 
advantage of fiexibility for arbitrarily unstructured meshes, with a compact stencil, and 
with the ability to easily accommodate arbitrary /ip -adapt ivity. For more details about 
DG scheme for convection dominated problems, we refer to the review paper [17]. The 
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DG method was later generalized to the local DG (LDG) method to solve the convection 
diffusion equation [16] and elliptic equations [I]. It is stable and locally conservative, 
which makes it particularly suitable to treat the Poisson equation. In our previous work 
[9], we proposed the first DG solver for (11.11) and showed some preliminary numerical 
calculations for one- and two-dimensional devices. In this paper, we will carefully formulate 
the DG-LDG scheme for the Boltzmann-Poisson system and perform extensive numerical 
studies to validate our calculation. 

This paper is organized as follows: in Section 2, we review the change of variables in 
[26| [5]. In Section 3, we study the DG-BTE solver for ID diodes. Section 4 is devoted to 
the discussion of the 2D double gate MOSFET DG solver. Conclusions and final remarks 
are presented in Section 5. Some technical implementations of the DG solver are collected 
in the Appendix. 

2 Change of variables 

In this section, we review the change of variables proposed in [26l H] . 

For the numerical treatment of the system (11. ip . (11.31) . it is convenient to introduce 
suitable dimensionless quantities and variables. We assume Tl = 300 K. Typical values for 
length, time and voltage are i^, = 10~^m, t^, = 10~^^s and V^, = 1 Volt, respectively. Thus, 
we define the dimensionless variables 

X t ^ V ^ ^ , E 

{x,y,z) = —, ^ = \7 ' ^' " Y 

with = 0.1 KC^ and 

In correspondence to [SB] and [S], we perform a coordinate transformation for k according 
to 

k = ^'^'^ ^bTl ^^i^^ axw) ^/i, a/1 - cos v?, v^l - /i^ sin ip^ , (2.6) 
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where the new independent variables are the dimensionless energy w = , the cosine of 

the polar angle fi and the azimuth angle ip with ax = ksTia. The main advantage of the 
generalized spherical coordinates fl2.6p is the easy treatment of the Dirac distribution in the 
kernel (11.51) of the collision term. In fact, this procedure enables us to transform the integral 
operator (11. 4p with the nonregular kernel S into an integral-difference operator, as shown in 
the following. 

We are interested in studying two-dimensional problems in real space but, of course, in 
the whole three-dimensional k-space. Therefore, it is useful to consider the new unknown 
function $ related to the electron distribution function via 

X, y, w, /i, (f) = s{w)f{t, X, k) , 

where 

s{w) = \/ w{l + aK'w){l + 2aKw), (2.7) 

is proportional to the Jacobian of the change of variables (12.61) and, apart from a dimensional 
constant factor, to the density of states. This allows us to write the free streaming operator 
of the dimensionless Boltzmann equation in a conservative form, which is appropriate for 
applying standard numerical schemes used for hyperbolic partial differential equations. Due 
to the symmetry of the problem and of the collision operator, we have 

$(t, X, y, w, /i, 2tt-(p) = $(t, X, y, w, /i, if) . (2.8) 

Straightforward but cumbersome calculations end in the following transport equation for $: 

^ + -^ig^<^) + ^(^2$) + 1^(^73$) + 1^(^74$) + ^(^?5$) = C($) . (2.9) 

The functions gi {i = 1, 2, .., 5) in the advection terms depend on the independent variables 
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w, n, if as well as on time and position via the electric field. They are given by 



2c, 



1 + 2aKW ' 
1 — n'^ -^^(1 + axw) cos ip 
1 + 2aKW 
a/w(1 + axw^ 
1 + 2aKW 



fiE^(t,x,y) + a/1 - /i2 cos ip Ey{t,x,y) 



Ck 



95{-) = Ck 



sin If 



1 - jj."^ Ex{t,x,y) - fi cos (pEy{t,x,y) 
Ey{t,x,y) 



with 



U 2kBTL ^ 
Cx = -r\ — and Ck 



t^^qE^, 



L V m* y/2m*kBTL 
The right hand side of (12. 9p is the integral-difference operator 

C($)(t, X, y, w, /i, ip) = s{w) jco dip' j dfi' ^{t, x, y, w, /i', p') 

dip' J djj! [c+$(t, x,y,w + 7, /i', p) + c_<l>(t, x,y,w - 7, fi', p')] 
- 2'k[cos{w) + c+s{w - 7) + c-s{w + 7)]$(t, x, y, w, n, p) , 



where 



(co,c+,c_)= ^.^ * ^/2m* kBTL{Ko,{ng + l)K,nqK) , 7- 



hp ksTL 

are dimensionless parameters. We remark that the 5 distributions in the kernel S have been 
eliminated which leads to the shifted arguments of $. The parameter 7 represents the jump 
constant corresponding to the quantum of energy hwp. We have also taken into account (12.81) 
in the integration with respect to p' . Since the energy variable uj is not negative, we must 
consider null $ and the function s, if the argument — 7 is negative. 
In terms of the new variables the electron density becomes 



7^(t^t, ^^y) — / f (j^*^^ ^^y^ k) d\^ 



\/2m*k^T£ 

h 



Pit,x,y) 



where 



p{t,x,y)= / dw dfi dp^{t,x,y,w,fi,p) 
Jo J-i Jo 



(2.10) 
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Hence, the dimensionless Poisson equation writes 



d f d f 

dx \ ^ dx J dy \^ dy 



Cp [p{t,x,y) -N'D{x,y)] 



with 



MD{x,y) 



No{Lx, i^y) and Cp 



y^2m*kBTL\ ilq 



Choosing the same values of the physical parameters as in [26], we obtain 



(2.11) 



Co 


^ 0.26531 




^ 0.16857 


Cp ~ 


1830349. 


c+ 


^ 0.50705 


Ck 


^ 0.32606 




10. 


c_ 


^ 0.04432 
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2.43723 


Ok - 


:i 0.01292 




= 11.7 











Moreover, the dimensional x-component of the velocity is given by 



dw dji I d(^gi{w,ii)^{t,x,y,w,fi,Lp) 
J-i Jo 



pit,x,y) 



the dimensional density by 



and the energy by 



1.0115 X lO^*^ X p{t,x,y) 



0.025849 X 



+00 l*\ 1*11 

dw dfx dip w ^{t, x,y,w, p, if) 
J~i Jo 



p{'t,x,y) 

In some simplified models, we consider our device in the x— direction by assuming that 
the doping profile, the potential and thus the force field are only x— dependent. By cylin- 
drical symmetry, the resulting distribution function does not depend on ip. In this case, the 
Boltzmann transport equation is reduced to 
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(2.12) 
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with 



gi{-> = 



1 + 2aKW 



( \ o Vwjl + aKw) 

gs{-) = -2c, ^^^^^^ ^^E{t,x) 

1 -/i^ 



and 



C($)(t, X, w, /i) = |co7r J dfi' ^(t,x,w, fi') 

+n J dn' [c+^{t,x,w + 'J, n') + C-^{t,x,w - 'J, n')] 

— 27i[cos{w) + c+s{w — 7) + C-s{w + 7)]$(t, X, w, /i) , 

In terms of the new variables the electron density becomes 

n(t^t,i^x) = f k) rfk = f '^'^^ J^^ Zl\ p(t,x), 

Jr' \ ^ J 

where 

/+00 i-i 
dw / dfi ^{t, x,w, fi) . (2.13) 

Hence, the dimensionless Poisson equation writes 

r 9^ j = f^^^' " -^^^"^^^ ^^-^^^ 

and 

E = -Cv^. 
ox 

3 DG-BTE solver for ID diodes simulation 

We begin with formulating the DG-BTE solver for ID diodes. These examples have been 
thoroughly studied and tested by WENO in [1]. 

The Boltzmann-Poisson system fl2.12p and (12.141) will be solved on the domain 

xe[0,L], we[0,wmax], fx E [-1,1], 



where L is the dimensionless length of the device and wmax is the maximum value of the 
energy, which is adjusted in the numerical experiments such that 

X, w, /i) ~ for w > Wmax and every t, x, ii. 



In f l2.12p . gi and are completely smooth in the variable w and yU, assuming E is given 
and smooth. However, (74 is singular for the energy w = 0, although it is compensated by 
the s{w) factor in the definition of $. 

The initial value of / is a locally Maxwellian distribution at the temperature Tx,, 

$(0,x,w;,/i) = s{w)ND{x)e-'^M 

with the numerical parameter Ai chosen so that the initial value for the density is equal to 
the doping N]j{x). 

We choose to perform our calculations on the following rectangular grid, 



ikm, 



(3.15) 



where 2 = 1,... A^^., /c = 1, . . . N^iui m = 1^ . . . N^, and 



Xj_|_i — Xi lb 



It is useful that we pick A*"^ to be even, so the function gi will assume a constant sign in each 

cell i^^fcm* 

The approximation space is thus defined as 



(3.16) 



where P^{Qikm) is the set of all polynomials of degree at most £ on Qikm- The DG formulation 
for the Boltzmann equation (12.121) would be: to find G V/, such that 



(^4$h {vf,),dn + F+-F-+F;^ 

ikm 

c{^h)vhdn. 



F~ + - F~ 



(3.17) 
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for any test function Vh G V/. In fl3.17p . 



'i', I 1 /"A* I 1 



= / ^i$t;+(Xi„i,w,/i)rftf;d/i, 



a; , 1 ru , i 



X. I J fl 1 

X. , 1 ri^ t 1 



X . 1 J II 1 



, 1 ^ UI 1 
, 1 ^UI, 1 



IX. 1 ^UI, 1 



where the upwind numerical fluxes $, $ are chosen according to the following rules, 

• The sign of gi only depends on /i, if /i^ > 0, $ = otherwise, $ = $+. 

• The sign of only depends on fiE{t, x), if ^mE{t, Xi) < 0, $ = otherwise, $ = 

• The sign of (74 only depends on E(t,x), if E{t,Xi) < 0, $ = $~; otherwise, $ = 

At the source and drain contacts, we implement the same boundary condition as proposed 
in [B] to realize neutral charges. In the {w, /x)-space, non boundary condition is necessary, 
since 

• ai w = 0, gs = 0. At w = wmax, ^ is machine zero. 

• At /X = ±1, g^ = 0, 

F^,F~,Fl^,F~ are always zero. This saves us the effort of constructing ghost elements in 
comparison with WENO. 
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The Poisson equation (12.141) is solved by the LDG method on a consistent grid of (13.151) 
in the x— direction. It involves rewriting the equation into the following form, 

(3.18) 



q -- 

d 

— {erq) = R{t,x 



where R(t, x) = Cp [p(t, x) — Afnix)] is a known function that can be computed at each time 
step once $ is solved from (13.171) . and the coefficient er here is a constant. The grid we use 



is li 



* 2 ^^2 



with i = 1, . . . , Nx- The approximation space is 



with P^{Ii) denoting the set of all polynomials of degree at most £ on Jj. The LDG scheme 
for (13.181) is given by: to find q^, ^ V"/, such that 



qhVhdx 



[ ^h{vh).dx - i'hv^, (x,+ i ) + ^hV^{x,i) = 0, 
Jii 

- (^rqh{Ph)xdx + e^f^pl{x^^i) -cirq,^pl{x^_i) = / R{t,x)phdx (3.19) 
Jii 2 2 Jii 

hold true for any Vh,Ph ^ W^. In the above formulation, the fiux is chosen as follows, 

"llh = "^h, = ^r-qt - \^h\, where [^!h] = ^+ - At x = L we need to fiip the fiux 

to = ^t,i — ^r-qh ~ [^/i] ^'^ adapt to the Dirichlet boundary conditions. Solving 

(I3.19p . we can obtain the numerical approximation of the electric potential and electric 

field Eh = —c^qh on each cell /j. 

To summarize, start with an initial condition for the DG-LDG algorithm advances 

from to t""*"^ in the following steps: 

Step 1 Compute ph{t,x) = n J^^dw J^-^^dfi ^hit, x,w, p). 

Step 2 Use ph{t,x) to solve from (I3.19P the electric field, and compute Qi, i = 1,3,4. 
Step 3 Solve ( 13.17P and get a method of line ODE for 

Step 4 Evolve this ODE by proper time stepping from to if partial time step is 
necessary, then repeat Step 1 to 3 as needed. 
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We want to remark that, unlike WENO, the DG formulation above has no restriction 
on the mesh size. In fact, nonuniform meshes would be more desirable in practice. For 
small semiconductor devices, it is common to have nonsmooth doping profiles and strong 
applied electric fields. The nonsmooth doping profile will create a distribution function / 
with high densities in some regions but low densities in other regions. Only a nonuniform 
grid may guarantee accurate results without using a large number of grid points. The strong 
electric fields give high energy to the charge particles. Hence the distribution function 
has, for some fixed points in the physical domain, a shape that is very different from the 
Maxwellian distribution (the equilibrium distribution in the absence of electric field), see 
for example Figures 13.111 to I3.14[ Moreover, taking into account the exponential decay of / 
for large value of |k|, a nonuniform grid can save us tremendous amount of computational 
time without sacrificing accuracy of the calculation. In the following simulations, we use a 
nonuniform mesh and refine locally near the junction of the channel and near /i = 1 where 
most of the interesting phenomena happen. 

We consider two test examples: Si —n — n^ diodes of a total length of 1 and 0.25/im, 
with 400 and 50nm channels located in the middle of the device, respectively. For the 400?7,m 
channel device, the dimensional doping is given by A'^^ = 5 x 10^'' cm~^ in the n"*" region and 
No = 2 X lO^^cm"'^ in the n~ region. For the 50nm channel device, the dimensional doping 
is given by Nij = 5 x lO^^cm"'^ in the region and = 1 x lO^^cm"^ in the n~ region. 
Both examples were computed by WENO in [5]. 

In our simulation, we use piecewise linear polynomials, i.e. i = 1, and second-order 
Runge-Kutta time discretization. The doping A'^^ is smoothened in the following way near 
the channel junctions to obtain non-oscillatory solutions. Suppose No = N^ in the n"*" 
region. No = Nj^ in the n" region and the length of the transition region is 2 cells, then the 
smoothened function is (A''^ — A''^)(l — y^)^ + N^, where y = {x — Xq + Ax) / (2 Ax + 10"^*^) 
is the coordinate transformation that makes the transition region [xq — Ax, xq + Ax) varies 
from to 1 in y. 
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The nonuniform mesh we use for 400nm channels is defined as follows. In the x-direction, 
if X < 0.2 or X > 0.4, Ax = 0.01. In the region 0.2 < x < 0.4, Ax = 0.005. Thus, the total 
number of cells in x direction is 120. In the ty-direction, we use 60 uniform cells. In the 
/i-direction, we use 24 cells, 12 in the region /i < 0.7, 12 in yU > 0.7. Thus, the grid consists 
of 120 X 60 X 24 cells, compared to the WENO grid of 180 x 60 x 24 uniform cells. 

We plot the evolution of density, mean velocity, energy and momentum in Figure 13.11 
The solution has already stabilized at t = 5.0 from the momentum plots. The macroscopic 
quantities at steady state are plotted in Figure [3^ The results are compared with the WENO 
calculation. They agree with each other in general, with DG offering more resolution and a 
higher peak in energy near the junctions. Figures [3.31 and 13^ show comparisons for the pdf 
at transient and steady state. We plot at different position of the device, namely, the left, 
center and right of the channel. We notice a larger value of pdf especially at the center of 
the channel, where the pdf is no longer Maxwellian. Moreover, at t = 0.5, xq = 0.5, the pdf 
shows a double hump structure, which is not captured by the WENO solver. All of these 
advantages come from the fact that we are refining more near /i = 1. To have a better idea 
of the shape of the pdf we plot $(t = 5.0, x = 0.5) in the cartesian coordinates in Figure 
13.91 The coordinate VI in the plot is the momentum parallel to the force field ki, V2 is 
the modulus of the orthogonal component. The peak is captured very sharply compared to 
WENO. 

The nonuniform mesh we use for 50nm channels is defined as follows. In the x-direction, 
near the junctions, in 0.09 < x < 0.11 and 0.14 < x < 0.16 , Ax = 0.001; in center of 
the channel 0.11 < x < 0.14, Ax = 0.005; at everywhere else. Ax = 0.01. Thus, the total 
number of cells in x direction is 64. In the w-direction, we use 60 uniform cells. In the 
/i-direction, we use 20 cells, 10 in the region /i < 0.7, 10 in /i > 0.7. Thus, the grid consists 
of 64 X 60 X 20 cells, compared to the WENO calculation of 150 x 120 x 16 uniform cells. The 
evolution and steady state plots are listed in Figures [3.51 to [3^ The conclusions are similar 
with 400nm, that we obtain better resolutions near the channel junctions and the peak for 
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Figure 3.1: Time evolution of macroscopic quantities using DG method for 400nm channel 
at V|gjg^g = 1.0. Top left: density in cm~^ ; top right: mean velocity in cm/s; bottom left: 
energy in eV; bottom right: momentum in cm~^ s~^. 

pdf is much higher. Figure 13.101 plots ^{t = 5.0, x = 0.125) in the cartesian coordinates. 
The peak is twice the height of WENO and is very sharp. Figure 13.111 to 13.141 plot the 
pdf near x = 0.15, the drain junction. We obtain distributions far away from statistical 
equilibrium, that reflects the lack of suitability of the classical hydrodynamical models for 
the drain region of a small gated device under even moderate voltage bias. 

We also compare the results from DG-BTE solver with those obtained from DSMC 
simulations, see Figures 13.151 13.161 The two simulations show good agreement except for 
energy plots near the boundaries. The modeling of the contact boundaries is not simple. 
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since it requires to know the distribution function of entering particles. The best way to solve 
this problem is the inclusion of a transport kinetic equation for the dynamic of the electron 
at the metal junctions; of course this is not realistic due to the complexity of this new 
kinetic equation, where the importance of electron-electron interaction requires a nonlinear 
coUisional operator, similar to the classical one of the Boltzmann equation for a rarefied 
perfect gas. Then, the simplest reasonable rule consists in assuming that the distribution 
function near, but outside, the device is proportional to a Maxwellian (or shifted Maxwellian) 
equilibrium distribution function, or to the distribution function near, but inside, the device 
boundaries. When there are strong electric fields also near the boundaries, the first choice is 
not reasonable, since, as we show in this paper, the distribution function is very far away to 
a Maxwellian distribution function. Therefore, the second choice is better than the first. We 
remark that both choices are simple but only low level approximation of the true physical 
phenomena; so many criticisms are known in the literature. We assume, as usual, that 
charge neutrality holds at the contact; so, the particle density near the contact boundaries 
coincides with the doping density. This law is used in all of the DSMC, WENO and DG 
simulations. Nevertheless, since we must approximate this constraint in different way, i.e. 
at molecular level for DSMC, introducing suitable ghost points for WENO scheme or giving 
appropriate values of $ at boundaries in DG simulations, we cannot have a unique exact 
boundary condition in the computational experiments. Now, it is obvious that this difference 
in the boundary treatment has an influence for the solutions at the stationary regime. 
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0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 

X X 

Figure 3.2: Comparison of macroscopic quantities using DG (symbols) and WENO (solid 
line) for 400nm channel at t = 5.0, Vj^jg^g = 1-0. Top left: density in cm~^; top right: mean 
velocity in cm/s; middle left: energy in eV; middle right: electric field in kV/cm; bottom 
left: potential in V; bottom right: momentum in cm~^s~^. Solution has reached steady 
state. 
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Figure 3.3: Comparison of the snapshot for ^{xo,w,fi) using DG (left) and WENO (right) 
for 400nm channel at t = 0.5, Vj^jg^g = 1.0. Top: Xq = 0.3; middle: Xq = 0.5; bottom: 
Xq = 0.7. Solution has not yet reached steady state. 
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Figure 3.4: Comparison of the snapshot for ^{xo,w,fi) using DG (left) and WENO (right) 
for 400nm channel at t = 5.0, Vj^jg^g = 1-0. Top: Xq = 0.3; middle: Xq = 0.5; bottom: 
Xq = 0.7. Solution has reached steady state. 



19 




Figure 3.5: Time evolution of macroscopic quantities using DG method for 50nm channel 
at Vj^ig^g = 1.0. Top left: density in cm~^; top right: mean velocity in cmjs\ bottom left: 
energy in eV\ bottom right: momentum in cm~^ . 
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Figure 3.6: Comparison of macroscopic quantities using DG (symbols) and WENO (solid 
line) for 50nm channel at t = 3.0, V^j^g^^ = 1.0. Top left: density in cm~^\ top right: mean 
velocity in cmjs\ middle left: energy in eV\ middle right: electric field in kV/cm] bottom 
left: potential in V\ bottom right: momentum in cm~^s~^. Solution has reached steady 
state. 
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Figure 3.7: Comparison of the snapshot for ^{xo,w,fi) using DG (left) and WENO (right) 
for 50nm channel at t = 0.5, V^jQig^g = 1.0. Top: Xq = 0.1 ; middle: Xq = 0.125; bottom: 
Xq = 0.15. Solution has not yet reached steady state. 
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Figure 3.8: Comparison of the snapshot for ^{xo,w,fi) using DG (left) and WENO (r 
solution for 50nm channel at t = 3.0, Vj^jg^g = 1.0. Top: Xq = 0.1; middle: Xq = 
bottom: xq = 0.15. Solution has reached steady state. 
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Figure 3.9: PDF for 400nm channel at t = 5.0, x = 0.5. 




Figure 3.10: PDF for 50nm channel at t = 3.0, x = 0.125. 
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;ure 3.11: PDF for 50nm channel at t = 3.0, x = 0.149. 






26 




Figure 3.15: Comparison of macroscopic quantities using DG (dashed line) and DSMC (solid 
line) for 400nm channel at t = 5.0, Vj^jg^g = 1-0. Top left: density in cm~^; top right: mean 
velocity in cm/s; middle left: energy in eV; middle right: electric field in kV/cm; bottom 
left: potential in V; bottom right: momentum in cm~^s~^. Solution has reached steady 
state. 
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Figure 3.16: Comparison of macroscopic quantities using DG (dashed line) and DSMC (solid 
line) for 50nm channel at t = 3.0, V^j^jg^g = 1-0. Top left: density in cm~^; top right: mean 
velocity in cm/s; middle left: energy in eV; middle right: electric field in kV/cm; bottom 
left: potential in V; bottom right: momentum in cm~^s~^. Solution has reached steady 
state. 
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4 DG-BTE solver for 2D double gate MOSFET simu- 
lation 

In this section, we consider a 2D double gate MOSFET device. In order to solver the 2D 
Boltzmann-Poisson system, we choose to implement a simple rectangular grid and let 



ijkmn 



X 



n+i 



where i = 1, . . . N^, j = 1, . . . Ny, A; = 1, . . . A''^, m = 1, . . . A''„, n = 1, . . . AL, and 



^i±\ — ± 2 



_ A/i^ 



The approximation space is defined as 



Av3„ 



(4.20) 



Here, P^{Qijkmn) is the set of all polynomials of degree at most £ on ^lijkmn- The DG 
formulation for the Boltzmann equation (12. 9p would be: to find $/i G V/, such that 



/ 

Jn 



gi^h {vh)x dQ 



ijkmn 



g2^h {Vh)ydVt 



ijkmn 



gz^h{,Vh)yjdVL- / gi^h{,Vh)^,dVL- / g^^h{,Vh)^dVL (4.21) 

^ijkmn *^ ^ijfcmn *^ ^ijkmn 



+ - P- + - F~ + - F~ + - F~ + - F~ 

' X X ' y y ' w w ' ifi ifi 



c{<^>^,)vf,dn. 



£~2 j ^ i. . 



ijkm.n 



for any test function Vh G ¥[. In (I4.2ip . 



i pfc+i r/^^+i /■'/'„+i 



^o-^l "'/^m-l '''^n-l 

'"^ "3 



Vj + ^ f^k+i /■^n+l 



fi-i ^ Vh (^i+i ^y^w.ix, ip)dy dw dfi dip, 
gi^v^{Xi_i,y,w,iJ,,ip)dydw dfi dip, 



■^^k-'^ "^/^m-l "^"^n-l 



- 1 1 /""'i I 1 /"Mil rf , 1 



X. I J W X J 1 ^l^ 1 



(72 "i'/i (a:, Vj^iyW, yU, (y9)(ix dwdfi dip, 
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r^+h r^+i r^+i - 

Fy = / / / g2^vj^{x,yj_i,w,fi,ip)dxdwdij.dip, 

J X ^ J W, 1 J U 1 J ifi 1 



'X. 1 ^«), 1 J II 1 1 



s; ,! rv , 1 /"M ii /"V 1 1 



""yj-l "'M™-! ""/'n-l 

^' ■ I 1 rv ■ , 1 /"A* I 1 , 1 



fl'3 (a;, 1 , fx, ip)dx dy dfi dip, 



^j-^ '^'"-2 '^""h 



f^i+i fVj+i f^k+^ fV„+i ^ 

= / / / / g4^Vi^{x,y,w,fi^^,(p)dxdydwdip, 
-^^i.i -^Vj^i J^k-i Jv^^i 

r^+h r^+h r^+i r-+i — 

Ff, = / / / g4^v^{x,y,w,iJ.^_i,(p)dxdydwdip, 

J X . 1 «/« . 1 J W, 1 J iO 1 



' X. I <'W. 1 1 >''P 1 



F^ = / / / g5^Vi^{x,y,w,fi,(p^^)dxdydwdfi, 

/•'"fc+i 

F^ = / / / 95^v^{x,y,w,fi,ip^_i)dxdydwdfi. 

Jx. 1 1 ^lii, 1 ^ (t 1 

1-1. f'-i "i-^ 

where the upwind numerical fluxes $, $, (74 $, $ are defined in the following way, 

• The sign of gi only depends on /x, if /x^ > 0, then $ = otherwise, $ = 

• The sign of g2 only depends on cose/?, if cos (fn > 0, then $ = $~; otherwise, $ = 
Note that in our simulation, is always even. 



For (73 $, we let 



^3 $ = -2c, 



k- 



nE^{t,x,y)^ + >/T^^ cos LP Ey{t,x,y)^ 



1 + 2aKW 

If firnEx{t,Xi,yj) < 0, then $ = otherwise, <l> = $+. 

If {cos (pn)Ey{t,Xi,yj) < 0, then $ = otherwise, $ = $+. 

For $, we let 

/1 r72 

a/1 - fj? Ex{t,x,y)^ - fj, cos if Ey{t,x,y)^ 
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If Exit,Xi,yj) < 0, then $ = otherwise, $ = $ 



If fimCOs{ipn)Ey(t,Xi,yj) > 0, then $ = $ ; otherwise, $ = $ 



• The sign of only depends on Ey(t, x, y), if Ey(t, Xi, yj) > 0, then $ = $ ; otherwise. 



The schematic plot of the double gate MOSFET device is given in Figure I4.17[ The 
shadowed region denotes the oxide-silicon region, whereas the rest is the silicon region. 
Since the problem is symmetric about the x-axis, we will only need to compute for y > 0. 
At the source and drain contacts, we implement the same boundary condition as proposed 
in [6] to realize neutral charges. A buffer layer of ghost points of z = and i = + 1 is 
used to make 



At the top and bottom of the computational domain (the silicon region), we impose the 
classical elastic specular boundary reflection. 

In the {w, fj,, v9)-space, no boundary condition is necessary, the reason is similar as in ID, 

• at w = 0, gs = 0. At w = wmax, ^ is machine zero; 

• at /i = ±1, = 0; 

• at = 0, TT, 5(5 = 0, 

so at the w,fi,(f boundary, the numerical flux vanishes, hence no ghost point is necessary. 

For the Poisson equation, \1/ = 0.52354 at source, \1/ = 1.5235 at drain and \1/ = 1.06 
at gate. For the rest of boundaries, we impose homogeneous Neumann boundary condition, 
i.e., 1^ = 0. The relative dielectric constant in the oxide-silicon region is = 3.9, in the 
silicon region is = 11.7. 




and 
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Figure 4.17: Schematic representation of a 2D double gate MOSFET device 



The Poisson equation (12. lip is solved by the LDG method. It involves rewriting the 
equation into the following form, 



^ dx ' 



dy 



^ (erg) + TT- (CrS) = R{t, X, y) 

ox oy 



(4.22) 



where R{t, x, y) = Cp [p(t, x, y) — Md^x, y)] is a known function that can be computed at each 
time step once $ is solved from (14.211) . and the coefficient depends on x,y. The Poisson 



* 2 



system is only on the (a;, y) domain. Hence, we use the grid lij 
with 2 = 1,..., j = 1, . . . , Ny + My, that includes the oxide-silicon region and is consistent 
with the five-dimensional rectangular grid for the Boltzmann equation in the silicon region. 
The approximation space is defined as 



(4.23) 



Here P^{Iij) denotes the set of all polynomials of degree at most i on Jjj. The LDG scheme 
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for (14.220 is: to find g^, Sh, e V/, such that 

/ qhVhdxdy+ / <I/h{vh)xdxdy - {x^j^i ,y)dy + / '^hvti^i-^^y)dy = 0, 

I ShWhdxdy+ I '^h{wh)ydxdy - I ^ ^ hW^{x,yj^i)dx + j ^ ^hW'^{x,yj_i)dx = 0, 



erqh{ph)xdxdy + I e^qf^p^{x.-^^,y)dy - I erqhVt{Xi_^,y)dy 
erSh{ph)ydxdy + / e^s^p^ (x, y^-^i )da; - / erSf,p+{x,yj_i)dx 

J X . I J X. I 

X, y)phdxdy (4.24) 

hold true for any Vh,Wh,Ph ^ W^^. In the above formulation, we choose the flux as follows, 
in the x-direction, we use = \E^^, e^i^ = erq^ — [^h\- In the ^/-direction, we use = "^h, 
e^Sh = ^r^h ~ [^h]- Near the drain, we are given Dirichlet boundary condition, so we need to 
flip the flux in X— direction: let ^h(Xj_,_i , y) = '^'l{x^^i,y) and e^^(Xj_,_i , y) = erq^{x,^^i,y) — 
[\l//i](Xj_,_i , ?/), if the point {x-^i,y) is at the drain. For the gate, we need to flip the flux in 
y-direction: let ^hix,yJ_^_l) = ^'^[(x, y^+i ) and e;:s;,(x, y^+i) = e^s;;(x, y^+i) - [^'/,](x, y^+i), 
if the point (x, Z/j+i ) is at the gate. For the bottom, we need to use the Neumann condition, 
and flip the flux in y-direction, i.e., = "^t^ = ^rSj^- This scheme described above will 
enforce the continuity of \& and across the interface of silicon and oxide-silicon interface. 
The solution of (14.241) gives us approximations to both the potential "ifh and the electric field 
iEx)h = —Cvqh, {,Ey)h = —CySh- 

To summarize, start with an initial condition for the DG-LDG algorithm for the 2D 
double gate MOSFET advances from to t""*"^ in the following steps: 

Step 1 Compute p/i(t,x,?/) = jQ°^dw J^^djj, J^dip ^h{t, x,y,w, i^, cp). 

Step 2 Use ph{t,x,y) to solve from (14.241) the electric field {Er^)h and {Ey)h, and compute 
gi,i = l,...,5. 

Step 3 Solve (14.211) and get a method of hue ODE for ^h- 
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Step 4 Evolve this ODE by proper time stepping from to if partial time step is 
necessary, then repeat Step 1 to 3 as needed. 

All numerical results are obtained with a piecewise linear approximation space and first 
order Euler time stepping. Apparently the collision term makes the Euler forward time 
stepping stable. We use a 24 x 14 grid in space, 120 points in w, 8 points in fi and 6 points 
in (f. In Figures 14.181 and 14.191 we show the results of the macroscopic quantities. We also 
show the pdf at six different locations in the device in Figure I4.20[ These pdf's have been 
computed by averaging the values of over if. In Figure 14. 2H we present the cartesian 
plot for pdf at {x,y) = (0.125,0.12), where a very non-equilibrium pdf is observed. 




Figure 4.18: Macroscopic quantities of double gate MOSFET device at t = 0.5. Top left: 
density in cm~^] top right: energy in eV; bottom left: x-component of velocity in cm/s; 
bottom right: y-component of velocity in cm/ s. Solution reached steady state. 
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Figure 4.19: Macroscopic quantities of double gate MOSFET device at t = 0.5. Top left: 
x-component of electric field in kV/cm; top right: y-component of electric field in kV/cm; 
bottom: electric potential in V. Solution has reached steady state. 
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Figure 4.20: PDF of double gate MOSFET device at t = 0.5. Top left: at (0.025, 0.12); 
top right: at at (0.075, 0.12); middle left: at (0.125, 0.12) ; middle right: at (0.1375, 0.06); 
bottom left: at (0.09375, 0.10) ; bottom right: at (0.09375, 0.) . Solution reached steady 
state. 
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Figure 4.21: PDF for 2D double gate MOSFET at t = 0.5, {x,y) = (0.9375,0.10). 
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5 Conclusions and final remarks 



We have developed a DG scheme for BTEs of type (11.11) . which takes into account optical- 
phonon interactions that become dominant under strong energetic conditions. We used the 
coordinate transformation proposed in [261 S] ^ind changed the colhsion into an integral- 
difference operator by using energy band as one of the variables. The Poisson equation is 
treated by LDG on a mesh that is consistent with the mesh of the DG-BTE scheme. The 
results are compared to those obtained from a high order WENO scheme simulation. By a 
local refinement in mesh, we were able to capture the subtle kinetic effects including very 
non-equilibrium distributions without a great increase of memory allocation and CPU time. 
The advantage of the DG scheme lies in its potential for implementation on unstructured 
meshes and for full /ij>-adaptivity. The simple communication pattern of the DG method 
also makes it a good candidate for the domain decomposition method for the coupled kinetic 
and macroscopic models. 
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A Appendix 

In this appendix, we collect some technical details for the implementation of the 2D DG-BTE 
solver. The discussion for ID solver is similar and omitted here. 

A.l The basis of the finite dimensional function space. 

In every cell ^ijkmn, we use piecewise linear polynomials and assume 

2(x — X '] '2{v — V ') 

, X, y, W, /i, ip) — ^jjfcmn(^) ~l" -^ijkirm(t) a ~^ ^ijkmnit) 



Ax, '^^'""^ ' Ay 



AWk Afirn Ay?„ 



It will be useful to note that 
^hit,x,y,w,fi,(p) = ^ 

k=l 



-Lijkmny'^) i ^ijkmny'^ ) ^ ijkmny'^ ) A?/ 

Xk{w) 



Awfe A/im Av9„ 

for every (x, ?/, w, /i, G flij^mn- Here, Xfc('W^) is the characteristic function in the interval 



fc=i 



A. 2 Treatment of the colUsion operator 

The gain term of the coUisional operator is 

G{^h){t,x,y,w) = s{w) \co dip' djj' ^h{t,x,y,w, fi'^ip') 



(. JO J-l 

fir fl 

+ I dip' d^' [c+^h{t,x,y,w + -f,fi',ip') + c-^h{t,x,y,w - -f,fi',ip')]\ . (A.2) 



J-l 



Now, we define 



{vh)^n{x,y,w) := / dip diJ,Vh{x,y,w, n,ip) 

' ip _i M _ 1 
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and, for a = —7, 0, 7, we have 



ijkmn 



dx dy dw dfi dip 



dip' / dfi' 



J-i 

EE 



/•TT f i 

s{w) I dip' I dfi' ^{t, x,y,w + a, fi' , ip' 

dy I dw s{w) x,y,w + a,fi ,ip) {vh)mn{x, y, w) 

Jw 1 

s{w) ^h{t,x,y,w + a, fi',ip') {vK)nin{x,y,w) dx dy dw dfi' dip' . 



dx 



m' = l n' = l ^'■ijkm'n' 

Now we discuss the following integral for different test function Vh, 



1=1 Vh{x,y,w,fi,ip) 



For Vhix,y,w,fi,ip) = 1, 



s{w) / dip' dfi' ^h{t,x,y,w + a, fi',ip' 
Jo J-i 



dx dy dw dfi dip . 



7V„ 

/ = ^ ^ ^ Afim' AiPn- 



m'=ln'=lk'=l 



+ Wijk'm'n'{t) / S{w) 



w._ 1 



Tijk'm'n' (^) 
2{w + a — Wk 

Awk' 



Xk'{w + a) dw 



Xk'{w + a) dw 



Axi Ayj Afim Av^n • 



2{x-Xi) 
Axi ' 



For Vh{x,y,w,fi,ip) 

I = -AXi Ayj Afim Aipn XI XI X] ^Vk'm'n'{t) 






X 



w. 1 



m'=l n'=l fe'=l 

Xk'{w + a)dw . 



For ^/^(x, ?/, w, fi, ip) = 



Ny__ 



/ = -AXi Ayj Afim Aipn ^Vk'm'n'{t) 



m'=l n'=l k'=l 



X / s{w) Xk' {w + a)dw . 
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For Vh{x, y, w, /x, (p) = 



A^M 

/ = ^ ^ ^ A/i„/ Aipn' 
m'=l n'=l k'=l 



Tijk'm'n' (^) 



w 1 



2 w - Wfc 

Xk'{u] + a) dw 

Awk 



+ ^ijk'm'n' (t) 



1 



4(w + cr-Wfc/)(w- Wfc) 

Xk'{w + a) dw 

Awk' Awk 



Axi Ay. A^rn Av^n . 



• For Vh{x, y, w, /x, ip) = ^%£f^, / = 0. 

• For Vh{x, y, w, /x, (p) = ^^^"^"^ , / = 0. 
The lost term in the colhsion operator is 



2Tr[cos{w) + c^s{w - 7) + C-.s{w + 7)]$(t, x, y, w, /i, (p) . 



(A.3) 



Let 



z/(w) = 27i[cos{w) + c+s{w — 7) + c_s(w + 7)] 



then we need to evaluate numerically, 



z/(w) X, y, w, fi, ip) Vh{x, y, w, /x, ip) dx dy dw dji dtp . 



ijkmn 



For Vh{x,y,w,fi,(p) = 1, 



I' = Axi Ayj Afirn Av^n 



X 



Tijkmnit) I U{w) dw + Wijknin{t) / 'u{w) 



h 2(w-Wk) 



Awk 



dw 



For Vh{x, y, w, v?) = 



/' = - AXi A^/j A/im A(y9„ Xijkmnit) 



z/(u;) (iu" . 
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For Vh{x, y, w, /x, (p) = 



u{w) dw . 



For Vh{x, y, w, (p) = ^^^^ 



I' = Axi Ay,j Afim. Aifn 



Tijkmn{t) I y{w) 



W 1 



2{w - Wk) 
Awk 



dw 



dw 



For w, /i, = "^^^^ 



I' = -AXi Ayj Aflm Aipn Mijkmnit) 



z/(w) dw . 



For ?/, w, /i, = 



I' = ^ AXi At/j A;U„ Aipn Pijkmnit) 



w. 1 
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A. 3 Integrals related to the collisional operator 

We need to evaluate (some numerically) the following integrals 

s{w) Xk'{w + a) dw 
2{w + cr — Wk') 



S[W] 



w, 1 



Xk'{w + a) dw 



^y^) ~. Xk'[w + cr) dw 

Wf,_ 1 



^Wk 



w, 1 



4(ii; + cT-tf;fc/)(u;- Wfc) , 
s[w) Xk' [w + cr) dw 



Awk' Awk 



uiw) dw 



, , 2(w - Wk) 



w, 1 



u^w) 



2{w - Wk) 
Awk 



dw 



If we evaluate these integrals by means of numerical quadrature formulas, then it is appro- 
priate to eliminate the singularity of the function s{w) at w = by change of variables. 



s{w — 7) dw 



s{w + 7) dw 



s{w) dw = / a/ w{1 + aKw){l + 2aKw) dw 



axr^ (1 + 2aKr )2r dr , {w = r ) 



6—7 
a— 7 



s{w) dw 



a/I + ai^r2 (1 + 2oLKr^)2r'^ dr , («; = + 7) 



'a— 7 
b+7 



s{w) dw 



a+7 
\/5+7 

'a+7 



1 + axr"^ {I + 2aKr ) 2r dr. {w = r — ^) 
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A. 4 Integrals related to the free streaming operator 



We recall that 

9i{-) = 

92{-) = 

9si-) = - 

9^{-) = - 

95{-) = Ck 
Now, we define: 



1 + 2aKW 
a/1 — iJ?^Jw{l + aKw) cos 



2ci 



Ck 



1 + 2aKW 
a/w(1 + aKw) 
1 + 2aKW 

sin (y9 



nE^{t,x,y) + a/1 - /i^ cos Ey{t,x,y) 



1 - jj."^ E^{t,x,y) - jj, cos (pEy{t,x,y) 
Ey{t,x,y) . 



1 + 2aii'W 
We need to evaluate the integrals. 



Si{w) dw 



S2{w) dw 



a/u;(1 + axw) 
1 + 2aKW 



a/w(1 + axw) 



dw 



a/1 + OiKr ^ 2 



^ 1 + 2ai^r2 



1 



dw 



2 r rfr , 



dr , 



, , 2(w; - Wfc) v'w^(l + Oi^w;) 2(u; - Wk) , 

dw = / — dw 

Awk J a 1 + 2aKW Awk 

^'VTT^2{r^-Wk)^^,^^ 



1 + 2aKr'^ ^Wk 



. ^2{w-Wk) . f' 1 2{w-Wk) . 

S2[w) dw = — dw 

Awk Ja ^yw{l + aKw) AWk 

^ VTTa^ Awk 
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'2{w - WkY 


la 


Awk 



1 2 



dw 



a 1 + 2aKW 



2{w - Wk) 



Awk 



dw 



1 S2{w) 


'2{w - WkY 


J a 


Awk 



n 2 



dw 



^ 1 + 2aKr' 

b 



2(r2 - Wk) 



-\ 2 



1 



a a/w(1 + axw) 



Awk 
2{w - Wk) 



Awk 



2r'^ dr . 



dw 



2{r^-Wk) 
Awk 



1 2 



dr . 



A. 5 Initial condition 



Since measure{flijkmn) = Axi Ayj Awk Afim Aipn, we have 

/ <l>/,(0, X, y, w, /i, (p) dx dydw djj, d^p = measure{^ijkmn) ^ijfcmn(O) , 

{x Xj)^ djjidif = - measure{Vtijkmn) -^ijfcmn.(O) ; 

3 



/ 

Jo., 



I ^h{0,x,y,w,^,ip) — — 

— ^ dx dy dw dadio = - 
Ayj 3 

1 



z\y — 

^h{0,x,y,w,^,ip) — — 
o , Ay, 

— dx dy dw djidip = - measure 
Awk 3 



^meaiSUreiytijkmn) ^ijkmniS^) 1 



f ^ .2{w-n 
j ^h{0,x,y,w,n,ip) — - — 

— — dx dy dw d^ dip = -measureiVtijkmn 

Aflm 3 
2(<y9 — (fn) 1 

$/i(0, X, y, w, n, if) dxdydw dn dp = -measure 

Apn 3 



/ 



^h{0,x,y,w,^,p) — — 



'(S^ijkmn) ^^ijkmni^^) ; 

Mijkmni^) , 
'■(S^ijkmn) -PijkmniS^) ■ 



If, for each {x, y) and (/i, p) , 

{F{x, y) s{w) for w < Wj^ 
otherwise 



and 



F{x,y) = Fij (constant) W{x,y) G 



1 • X,:\ 1 

* 2 '^2 



X 
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then it is reasonable to assume 



Tijkmn (0) 



s{w)e dw . 



W 1 



WijkmniO) = 3 



Awk 



^i, 1 



s(w)e T aw, 



Mijkmn{0) — PijkmniO) — 



Recalling the definition f l2.10p of the dimensionless charge density, we have 

p(0, = 2 TT y) ^ / " s{w) e-^ dw , 

fc=i -^^k-^ 

which gives a relationship between F{x, y) and the initial charge density p. 
A. 6 Hydrodynamical variables 



If {x,y) e 



then 



N^, 



Ph{t,x,y) = ^J2Y1 



k=l m=l n=l 



2{x - Xi) 



Axi 

AWkAprn Aipr. 



We define 



Afu, ^V' 

Tijit) = X] Tijkmn{t) AWkAprn Aifn 

fc=l m=l n=l 

TV™ 

fc=l m.=l n=l 
7V„ A^M Afv' 
= ^^^^^^yijkmnit) AWkAHraAifn- 



k=l m=l n=l 



Therefore, for every [x, y) G 



^ / N / N 2(a; — Xj) - , , 2(y — u,) 



AaJi 



Ay, 



(A.4) 

(A.5) 
(A.6) 

(A.7) 



(Ai 



(A.9) 
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For every (x, y) G 



' 2 2 



the approximate momentum in x-direction is 



"y^^ A(/5„ Tijkmnit) + Qlw, km Wijkmn{t) + Ql^^km Mijkmn{t)\ 

fc=l m=l n=l 

A^J^ (71,^^ ^ijkmnit) 

_k=l m=l n=l 
' 

+ ^V^n 5'l,fcm yijkmn{t) 

k=l m=l n=l 



2 ^ 2/ ^2 
Axi 

Ay, 



(A.IO) 



and the approximate momentum in ^/-direction is 



Nu, 

^ ^ ^ ^ ^ ^ [92,kmn Tij/^mnii") ~l~ 92w,kmn ^^ijkmnii^ ~l~ 92^,kmn -^{jkmni^^") ~l~ 92ip,kmn Pijkmri{j^)\ 
k=l m=l n=l 

+ 92,kmn Xijkmn{t) 



k=l m=l n=l 
■ iV„ A^M 



2 fx — X,- 



^ ^ ^ ^ ^ ^ 92,kmn "^ijkmnit^ 



k=l m=l n=l 



Axi 



(A.ii: 



where 



91, km — 
9lw,km ~ 
9llJ,,km ~ 

92, kmn ~ 
92w,kmn 
92fj.,kmn 
92ip,kmn 



w, I 1 ru ,1 



9i{w, fi) dw dfi 



W, I J U 1 



"^i I 1 /"M I 1 



w, 1 J 1 

'"ill rfj- I 1 



W. 1 J II 1 

I 1 /"M I 1 /"V I 1 



2M;-U7fc 

5-1 [w, M T dw dfi 

Awk 

9i [w, fx) — dw dfx 

92{w, ^,ip) dw djji dip 



W, I J U I J Ifi 1 



w, \ J a \ 'J ^ 1 



2{w -Wk) 

92[w, /i, dw djji dip 

Awk 

92{w, /i, p) ^""^ (iu7 dfi dip 

Afim 



-,+ 1 /-M^+i f^„^i 2{ip-ip^) 

92\w, /i, dw djjL dip . 

Aipn 



w I -J fl x J ip 1 
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Analogously, the energy multiplied by the charge density ph(t,x,y) is 



Af„ 

k=l m=l n=l 

■ Nni -A^M '^v 



Wk AWk Tijkmn{t) + ^^^^^ Wijkmn{t) 



^ ^ ^ AV9„ A/im Wfc AWk Xijkmn{t) 
k=l m=l n=l 

Nu, ^-P 

^ ^ ^ A(y9„ A/im Wfc Az/Jfc Yijkmn{t) 



k=l m=l n=l 



6 

2(a: - Xj) 
Axi 

2(y - y,) 

A%- ' 



(A.12) 



for every (x, y) G 
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